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ABSTRACT 

In this • '■diminary study involving advanced CFD codes, an incremental for- 
mulae known as the “delta” or “correction” form, is presented for solving 

sparse systems of linear equations which are associated with aerody- 
ensitivity analysis. For typical problems in 2D, a direct solution method can 
be applied to these linear equations in either the standard or the incremental form, 
in which case the two are equivalent. Iterative methods appear to be needed for 
future 3D applications, however, because direct solver methods require much more 
computer memory than is currently available. Iterative methods for solving these 
equations in the standard form result in certain difficulties, such as ill-conditioning 
of the coefficient matrix, which can be overcome when these equations are cast in 
the incremental form; these and other benefits are discussed herein. The method- 
ology is successfully implemented and tested in 2D using an upwind, cell-centered, 
finite volume formulation applied to the thin-layer Navier-Stokes equations. Re- 
sults are presented for two laminar sample problems : 1) transonic flow through a 
double-throat nozzle, and 2) flow over an isolated airfoil. 
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1.0 Introduction 


For many complex flow fields of interest in practical engineering problems, accurate detailed 
analyses are now possible using supercomputers and advanced software; these codes have been 
developed in recent years through an intensive research effort focused in the discipline now 
known as Computational Fluid Dynamics (CFD). For these advanced CFD codes to become 
more useful as practical design tools, additional software is needed which will efficiently provide 
accurate aerodynamic sensitivity derivatives which are consistent with the discrete flow solutions 
of the particular CFD code of choice. The theme of this report is the ongoing development of 
a methodology for calculating these derivatives. 

A sensitivity derivative is defined as the derivative of a system response of interest (e.g., 
the lift or drag of an airfoil) with respect to an independent design variable of interest (e.g., a 
parameter which controls the shape of an airfoil). In a typical design environment, a very large 
number of analyses are often made in determining the “best” design. An efficient method for 
calculating accurate sensitivity derivatives can be applied in several different ways to significantly 
reduce the number and/or computational cost of these multiple analyses. This could be critical 
for the integration of advanced CFD codes into a systematic design methodology, where the 
computational cost of a single flow analysis can be extremely high, particularly in 3D. 

One method of a very general yet conceptually simple nature for computing aerodynamic 
sensitivity derivatives is the method of “brute force” finite differences. With this method, 
assuming forward finite difference approximations are used, the CFD flow analysis code is used 
to generate one converged flow solution for a slightly perturbed value of each design variable for 
which sensitivity derivatives are required. The principal drawback of this method is clearly that 
of computational cost, since the number of flow analyses required in a typical design problem can 
be extremely (i.e., prohibitively) large, particularly when the number of design variables is large. 

As a typically less costly alternative to the finite difference approach, aerodynamic sensitivity 
derivatives can (in principle) be computed by direct differentiation of the governing equations 
which control the fluid flow. If the continuous governing equations are differentiated prior to 
their numerical discretization, the method is known as the “continuum” approach. In contrast, 
if the resulting algebraic equations which model the governing equations are differentiated 
following their discretization, the method is known as the “discrete” approach. In developing 
efficient methods for computing these sensitivity derivatives and their subsequent application to 
aerodynamic design problems, researchers have been and remain active; Refs. 1 through 24 are 
a representative (but not exhaustive) sample of articles which are germane to the present effort. 
Reference 8 addresses the distinction between the aforementioned “continuum” and “discrete” 
approaches, and Ref. 24 is a concurrent study which addresses related issues of specific interest 
here. 

The present study represents an extension of the recent efforts of Refs. 13 through 23, where 
using the discrete approach, fundamental sensitivity equations are derived by direct differentiation 
of the system of nonlinear algebraic equations which model either the Euler or thin-layer Navier- 
Stokes (TLNS) equations for 2D steady flow. This differentiation results in very large systems 
of algebraic linear sensitivity equations which must be solved to obtain these derivatives of 
interest. In Refs. 13 through 23, the fundamental sensitivity equations are solved in what is 
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henceforth referred to herein as the “standard” (i.e., non-incremental) form. Furthermore, in these 
references, a direct solver method is applied to solve these equations; the single exception is Ref. 
23, where a hybrid direct solver/conventional iterative approach is « lopted for an isolated airfoil 
example problem. There are some important advantages in using a direct method when feasible; 
these are discussed in the references and also noted in a later section of this report. However, 
the most serious disadvantage of a direct solver method is the extremely large computer storage 
requirement, which for practical 3D problems appears to be well beyond the current capacity of 
modem supercorr outers; this capacity can even be exceeded in 2D on very fine grids. 

In an effort 10 circumvent the computer storage limitation for the direct methods, this 
preliminary study focuses on fundamental algorithm development for the efficient iterative 
solution of the aerodynamic sensitivity equations. That is, the print al motivation and objective 
is to develop a solid framework in 2D from which future extensions to 3D will be feasible. In 
gene; one of the most serious difficulties encountered in the development and/or application 
of iterative techniques is that of poor overall conditioning and lack of diagonal dominance 
in the coefficient matrix. Unfortunately, this is a very common occurrence in the coefficient 
matrices of interest here; the severity varies greatly and depends on many factors. This problem 
can manifest itself in either poor performance or even complete failure (i.e., divergence) of an 
iterative algorithm. 

A computationally useful property of the “incremental” form (also commonly known as the 
“delta” or “correction” form) can be effectively exploited to combat these problems of poor matrix 
conditioning. This property is that “approximations of convenience” can be introduced into the 
coefficient matrix of the equations, without affecting the final converged values of the sensitivity 
derivatives. The approximations must be “reasonable” enough that the resulting iterative strategy 
is convergent. In contrast, if any approximations are made to the coefficient matrix of the 
equations in the standard form, then the computed sensitivity derivatives cannot be consistently 
discrete; tha: is, they will not be the correct derivatives of the algebraic equations which are 
solved when generating the steady-state flow solution. In the implementation of the incremental 
formulation herein, a judiciously selected block-diagonally dominant matrix is introduced as 
an approximate rep ement for the original ill-conditioned left-hand side coefficient matrix. 
The positive impact which this can have on the development of iterative techniques for the 
aerodynamic sensitivity equations is discussed herein, and illustrated in the example problems. 
Additional benefits which might be derived from this flexible nature of the “delta” formulation 
are also discussed. 

The remainder of this report is organized as follows. The next section, presentation of theory, 
is further subdivided into five subsections which review and discuss: 1) the governing equations, 
2) the spatial discretization and implicit formulation, 3) the fundamental sensitivity equations 
in standard form, 4) basic linear equation solving in incremental form, and 5) incremental 
solution of the equations of sensitivity analysis, where some significant implications of this 
formulation compared to the standard form are noted. Following the presentation of theory 
section, computational results are presented, illustrating application of the methodology to two 
laminar viscous flow example problems: 1) transonic internal flow through a double-throat nozzle 
and 2) external flow over an isolated airfoil. The last section is a summary where conclusions 
are given. In an appendix, the direction of ongoing and future work is discussed, where sample 
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results are shown from the successful application of a spatially-split approximately factored 
strategy for efficiently solving the sensitivity equations in incremental form. 

2.0 Presentation of Theory 


2.1 Governing Equations 


The governing equations in this study are the 2D thin-layer Navier-Stokes (TLNS) equations; 
they are 


IdQ 
J dt 
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The vector, R(Q), is known as the residual, and is clearly null for steady flow. The elements of 
the vector, Q, are the conserved variables, where, p is density, u and v are velocity components 
in Cartesian coordinates, and eo is total energy (i.e., eo = e + " |~ v , where e is the specific 
internal energy of the fluid). The inviscid flux vectors, F(Q) and G(Q), are 

F(Q) = ^F(Q) + ^G(Q) 

v i (4) 

G(Q) = -j-F(Q) + ^G(Q) 


A transformation to generalized coordinates from Cartesian (x,y) coordinates has been 
made in Eq. (1), where £ x , £y , ?/x, Vy are “metric” terms, and J is the determinant of the Jacobian 
matrix of this transformation. The Cartesian flux vectors, F(Q) and G(Q), are 

T 

F(Q) = [/>U, pu 2 -(- p, puv, (peo + p)u] 
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G(Q) = [/9V, /)UV, pv 2 + p, (peo + p)v] 

The pressure, p, is evaluated using the ideal gas law 


P = (7 - 1) 
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and 7 is the ratio of specific heats, taken to be 1.4. The thin-layer viscous terms in generalized 
coordinates are 
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The molecular viscosity is given by n, Stokes’ hypothesis for the bulk viscosity (A = 
-2fi/3) has been used, a is the speed of sound, Pr is the Prandtl number (taken to be 0.72), 
and ReL is the Reynolds number. Nondimensionalization of Eq. (1) is with respect to 
and Uoo, the freestream density and velocity, respectively. The physical coordinates (x,y) are 
nondimensionalized by a reference length, L, and the viscosity is nondimensionalized by /Xqo, the 
molecular viscosity of the freestream. The nondimensional molecular viscosity can be computed 
using Sutherland’s law and a reference temperature. Too. the static temperature of the freestream. 
For additional simplicity here, however, the molecular viscosity is taken to be constant, equal 
to that of the freestream. 


2.2 Spatial Discretization and Implicit Formulation 

Computationally, the TLNS equations are solved here in their alternative integral conserva- 
tion law form using an upwind cell-centered finite volume formulation (see Refs. 25 through 
31), where the residual at each cell is evaluated as a balance of inviscid and viscous fluxes 
across cell interfaces. Upwind evaluation of the inviscid fluxes is accomplished by upwind 
interpolation of the field variables, Q, from the approximate cell centers to the cell interfaces, 
where the flux-vector splitting procedure of van Leer (Ref. 32) is employed. In this study, 
third-order accuracy is used for the inviscid flux balance in the streamwise (£) and in the normal 
0?) directions. The finite volume equivalent of second-order accurate central differences is used 
to approximate the thin-layer viscous terms. This results in a higher-order accurate algebraic 
approximate representation of the residual at each cell in the domain. When assembled globally 
including all cells and boundary condition relationships, this can be expressed as 

{R(Q'» = {0} (9) 

where {Q*} is called the “root” (i.e., the steady-state value of the field variables). Therefore, Eq. 
(9) represents a large coupled system of nonlinear algebraic equations; thus finding a steady-state 
solution to the TLNS equations has been replaced (approximately) by the problem of hading 
the root, {Q* }, of this set of algebraic equations. In Eq. (9) and henceforth, the notation, ‘{ }’, 
indicates a global column vector. 
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The TLNS equations are discretized in time using the Euler implicit method, followed by a 
Taylor’s series linearization of the discrete equations in time about the known time level. This 
results in a large system of linear algebraic equations at each time step, which is 

([jSt]-[y]) { * A<,, “ {R ‘ (Q)} (10) 

{Q"+>}={Q»} + {"AQ} (U) 

n = 1, 2, 3, ... 

Equations (10) and (11) represent the fundamental implicit formulation for integrating the 
TLNS equations in time to steady-state. In these equations, ‘n’ is the time iteration index, 
and { n AQ} is the incremental change in the field variables between the known (n th ) and the 
next (n ,h +l) time levels. The matrix, , is diagonal, and contains the time term. The 

large Jacobian matrix, ^10 1 ; i s sparse and has a banded structure, with nine diagonals, the 
individual elements of which are 4x4 block coefficient matrices. In addition to its use in Eq. 
(10) above, this important Jacobian matrix plays another central role in this study, which will 
be shown later. 

In principle, Eq. (10) can be repeatedly solved directly (using Eq. (11) to update the 
field variables), as the solution is advanced to steady-state; for very large time steps, the direct 
method represents Newton’s root finding procedure for nonlinear equations. The direct method 
however is not necessarily the most efficient procedure with respect to overall CPU time (Ref. 
33), and the large storage requirements of the method make it infeasible in 3D. Therefore, 
more commonly, an iterative algorithm is selected for use in the repeated solution of Eq. (10). 
Popular choices of these iterative algorithms include approximate factorization (AF) (Ref. 34), 
conventional relaxation algorithms (Refs. 29, 30), the strongly implicit procedure (SIP) (Ref. 
35), and preconditioned conjugate gradient methods (Refs. 36, 37), to name a few. 

It is noted that Eqs. (10) and (11) are an incremental formulation for solving the nonlinear 
problem of Eq. (9). If convergence is achieved, the steady-state solution, {Q*}, only depends 
on what is implemented in the discrete formulation of the residual vector on the right-hand side 
of Eq. (10). It is also implied that this solution is independent of any approximations which are 
made in the coefficient matrix of Eq. (10). The final solution is also independent of the initial 
guess, and all transient solutions which are generated prior to convergence. 

For typical advanced CFD flow codes which employ the implicit time integration formulation 
of Eqs. (10) and (11), the following approximations are often seen in the coefficient matrix of 
Eq. (10) (the list is a representative but not exhaustive one): 

1) A first-order accurate upwind spatial discretization of the implicit terms is used, even 
though a higher-order accurate spatial discretization, either upwind or perhaps even 
central “differences” (Ref. 29), is used on the right-hand side of the equation. 

2) A consistently linearized treatment of the boundary conditions in “delta” form is typically 
neglected. In particular, a fully consistent treatment of the implicit terms resulting from 
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the “periodic” boundary conditions of “C” and “O” meshes and also of the implicit terms 
across the zonal interfaces of multiblock grids is not used. 

3) Only approximate solutions of Eq. (10) are actually generated at each time step with 
the use of iterative methods, in order that each time step is efficiently completed. 

The preceding examples and many others not mentioned are “approximations of convenience” 
and are made on the left-hand side of Eq. (10) in order to influence the nature of the 
resulting algorithm which is to be used in finding the solution. These may be introduced for 
computational simplicity of implementation or overall efficiency, or both. This flexibility of 
the vita formulation, which allows approximations to be introduced into the left-hand side 
co; ,ent matrix without influencing the final solution, can also be exploited in the solution of 
the linear aerodynamic sensitivity equations, as will be seen in subsequent sections. 

23 Fundamental Sensitivity Equations In Standard Form 

ider the vector, /3, whose elements are independent variables, typically called the design 
v^ none, or all of these variables may be related to the geometric shape of the 

boundary ^urtace of the flow problem of interest. Computationally, the geometric shape of the 
domain is defined by the mesh upon which calculations are made; the complete vector of (x,y) 
coordinates which defines the mesh is represented symbolically herein as {X}. For a steady-state 
solution, the discrete residual vector given by Eq. (9) is expressed now in the following form 

{R(Q-(/3),X(^),«} = {0) (12) 

where the explicit dependence of the discrete residual on the computational mesh, {X}, as well 
as its explicit dependence (if any) on /3 has now been emphasized. Direct differentiation of Eq. 
(12) ■■ h respect to /3k, the element of 3, yields 

dRHdQ*] [dRHdX] fdR] 

dQ\\ d/3 k f~ LaxJ + l^kj (13) 

Term 1 Term 2 Term 3 

Equation (13 "an exact derivative of the discrete algebraic residual vector; this procedure is 
l ' n in Refs. _ and 4 as the quasi-analytical method. The Jacobian matrix, , of Eq. (13) 
cal to that found in the fundamental implicit formulation for numerical time integration 
vt-q. (10)) of the TLNS equations, except that is evaluated at steady-state, {Q*}. It is thus 
well understood. The solution vector, j j , is the sensitivity of the complete vector of field 

variables with respect to the k 1 * 1 design variable. The matrix, , is the Jacobian matrix of the 
discrete steady-state residual vector with respect to the complete vector of (x,y) grid coordinates; 
it is documented in detail in Ref. 17. The vector, j, of Term 2, contains what is referred to 
here as the grid sensitivity terms; these are the sensitivity derivatives with respect to /3k of each 
‘x’ and ‘y’ coordinate point of the entire computational mesh. The treatment of the terms of the 
grid sensitivity vector is given special consideration in Refs. 18, 23, 38, and 39. The vector. 
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accounts for derivatives resulting from explicit dependencies (if any) of the residual 
vector on 0k, and additional discussion concerning this is found in Ref. 21. In the event that 
fa is not a design parameter for the geometric shape, then the second term of Eq. (13) will 
be zero, since the vector, j j, is then null. If fa is a geometric shape design parameter, its 
effect on the residual (Eq. (12)) will usually be felt only through the grid, and the final term 
of Eq. (13) will generally be zero. 

It is strongly emphasized that all boundary condition relationships must be treated in a fully 
consistent manner, and included in Eq. (13) above. Proper boundary condition treatment should 
be included in the Jacobian matrices, and f|r , as well as in the vector, ||^|. If accurate 
results are to be obtained using the present methods, it is critical that this is not neglected here 
as it often is in the fundamental implicit time integration formulation (i.e., Eq. (10)). Detailed 
documentation on the consistent treatment of the boundary conditions and its importance in 
these equations is found in Refs. 21, 22, and 23. 

Note that Eq. (13) is a linear system of equations which in principle can be solved directly 
for the vector, {-^ }. Of course, the solution of Eq. (13) must be repeated for each element of 

0 (i.e., for each design variable) for which sensitivity derivatives are desired. However, a single 
LU factorization of the coefficient matrix can be repeatedly reused for multiple solutions (i.e., 
for multiple design variables) in the forward and backward substitution operations. The reuse 
of the LU factorization can represent a substantial savings in computational work, particularly 
when the linear system of Eq. (13) and/or the number of design variables of interest is large. 

The solution of Eq. (13) for the vector, is not the final goal; rather, the sensitivity 

derivatives of some specific system responses are sought (e.g., for an airfoil, the sensitivities 
of the lift, drag, and moment coefficients might be required). Consider therefore the j* system 
response of interest, Cj, which in general can be functionally dependent on the steady-state field 
variables, {Q* }, the grid, {X}, and also explicitly on the design variables, 0; that is 

Cj = Cj(Q*(/?),X (/?),/?) (14) 

The total rate of change of the j* system response, Cj, with respect to the k* design variable, 
0k, is then given by 

dCj _ /aCj\ T fSQll , jaCj\ T f ax i ac, 

d/3 k \dQj \dlhS \axf \d/3 k J dfa (15) 

Term 1 Term 2 Term 3 

where in Eq. (15), Term(s) 2 and/or 3 could be zero, depending on the particular system response 
(Cj) and design variable (0k) of concern. Solution of Eq. (13) therefore provides the vector, 

j, which is needed in Eq. (15). Furthermore, for geometric shape sensitivity derivatives, 

the grid sensitivity vector, |^-|,ofEq. (13) is reused, if needed, in Eq. (15). Specific ancillary 
sensitivity relationships of the type given by Eq. (15) which are used in the present study for 
computing sensitivity derivatives of aerodynamic force coefficients are presented in Ref. 23. 
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On the left-hand side of Eq. (15) above, the notation for a total derivative has been used, 
indicating that the total rate of change of Cj with respect to /? k is included in the expression, 
and to distinguish it from the partial derivative term (Term 3) on the right-hand side of the 
expression. However, it should be understood that this derivative is still a partial derivative 
in the sense that Cj is in general a function of multiple independent design variables. For 
consistency, this notation will continue to be used throughout. 

A closely related alternative procedure for computing sensitivity derivatives, known as the 
adjoint variable approach, is easily developed using the relationships presented thus far. This 
begins by combining Eqs. (13) and (15) to yield 

ACj = l^i\ T l 2011 , 3% \ aCj 

dft, \dQJ 1%/ \ax/ \%J 3ft 



The adjoint variable vector, { Aj }, is arbitrary at this point, since the inner product of { A j } is 
taken with the null vector, from Eq. (13). Thus there is no net change from Eq. (15) to Eq. 
(16), since the entire additional term on the right-hand side of Eq. (16) is zero, for any and all 
{Aj}. Rearranging, Eq. (16) becomes 



The necessity of evaluating the vector, using Eq. (13) is eliminated for all /3 k by 

selecting the vector, {Aj}, such that the coefficient of j^-j in Eq. (17) is null. That is, 
selection of {Aj} which satisfies 


implies 



(18) 


<9R 

.0Q. 


iT 


{A 


>»-{§} 


(19) 


Therefore, following the solution of Eq. (19) for this particular choice of the adjoint variable 
vector, {Aj}, the sensitivity derivatives of Cj with respect to all /? k are computed by 


dC; 

dA 



+ {Aj} T 


■0R‘ 

dX 




( 20 ) 
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Solution of the linear system of Eq. (19) for { Aj } is analogous to the solution of Eq. (13) 
for in that the respective coefficient matrices are transposes of each other. A particular 

solution, {Aj}, is valid only for a specific system response, Cj, and thus solution of Eq. (19) 
must be repeated for each different system response of interest If Eq. (19) is solved directly, 

however, multiple solutions require only a single LU factorization of the coefficient matrix, 

' dc ' 


which is repeatedly reused for an unlimited number of right-hand side vectors, | j (i.e., for 
an unlimited number of different system responses of interest). 

It is simple to verify from the preceding equations, and significant to note, that each solution, 
of Eq. (13) for a particular design variable can be used for an unlimited number of 
different system responses. In contrast, however, each solution, {Aj}, of Eq. (19) for a particular 
system response can be used for an unlimited number of different design variables. Therefore, 
in terms of computational work, if the number of system responses of interest is larger than 
the number of design variables, then sensitivity derivatives should be computed by solving Eq. 
(13). Otherwise greater computational efficiency is obtained using the adjoint variable method. 
Despite this difference which has been noted between these two closely related procedures, it is 
emphasized that the two methods are equivalent in the sense that they yield identical values for 
the sensitivity derivatives, if properly implemented computationally. 

The significance of the well-known difference in the computational efficiency of the two 
methods is mitigated greatly if a direct method is used to solve the linear systems (i.e., either 
Eq. (13) or Eq. (19)), because the LU factorization must only be done once for multiple right- 
hand side vectors. However, this distinction becomes very important if an iterative strategy is 
used to solve these Unear systems, particularly if the difference between the number of design 
variables and the number of system responses of interest is very large. This difference occurs, 
of course, because with iterative methods, the computational work required for solution of each 
linear system is approximately equal to the computational work required to solve the first one. 

Summarizing briefly, it has been shown that calculating aerodynamic sensitivity derivatives 
using the discrete direct differentiation method requires the solution of large linear systems of 
equations of the type given by a choice of either Eq. (13) or Eq. (19). Henceforth in this 
report, these two systems of linear equations are known as the aerodynamic sensitivity equations 
in standard form . Fundamental algorithm development for the iterative solution of one of these 
two linear systems is easily extended and applied to the other, since as noted previously, their 
respective coefficient matrices are transposes of each other. In the example problems for which 
sensitivity derivatives are calculated in a later section, actual implementation and testing of the 
methods proposed herein is accomplished using Eq. (13), although the adjoint variable method, 
Eq. (19), could also have been used. When the linear aerodynamic sensitivity equations are 
solved in standard form, it should be noted that no approximations can be introduced into any 
of the terms, without simultaneously introducing error into the resulting sensitivity derivatives. 
In this form, the framework to support the development of iterative methods is thus rigid and 
restrictive. 

As a consequence of the preceding discussion, for the higher-order accurate upwind spatial 
discretization which is selected herein for the flow analysis, a consistent higher-order accurate 
upwind spatial discretization including a fully consistent treatment of all boundary conditions is 
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required in the left-hand side coefficient matrix of the sensitivity equations (in standard form). 
Furthermore, there is no “time term” added here to enhance each element of the diagonal, as 
seen (in contrast) in the implicit time integration formulation of Eq. (10). Unfortunately, the 
resulting coefficient matrix in this case is not diagonally dominant (Ref. 29), and consequently 
the computational performance of traditional iterative methods for the sensitivity equations in 
standard form is expected to be poor, or even fail. If the present methods were applied using a 
popular “central difference” discretization of the inviscid terms in the flow solver, the diagonal 
dominance of the resulting sensitivity equations would become far worse. Therefore, it is this 
particular difficulty (i.e., the lack of a sufficiently strong diagonal) and how it can be overcome 
which is of principal concern in the development of the incremental form of the equations in 
the following sections. 


2.4 Basic Linear Equation Solving in Incremental Form 

Consider the linear system of algebraic equations in the general form 

[A]{Z*} + (B) = {0} (21) 

where {Z*} is the solution vector. In treating the problem of solving Eq. (21), in essence a “root 
finding” problem, application of Newton’s method (traditionally used in root finding for nonlinear 
equations) to the linear problem yields the basic two-step iterative incremental formulation 

Step 1 - [A]{ m AZ} = [A]{Z m } + {B} (22) 


Step 2 {Z m+1 }= {Z m } + { m AZ} 

m = 1,2,3,.... 


(23) 


where ‘m’ is an iteration index, and { m AZ} is the incremental change in the solution from the 
known (111 th ) to the next (m th +l) iteration level. An initial guess, j Z 1 j, is required to begin 
the procedure, which in the present study is taken everywhere as zero. If Newton’s method is 
applied strictly, the coefficient matrix [A] is equal to the matrix [A], and clearly the two-step 
iterative strategy of Eqs. (22) and (23) for the linear problem converges on the first iteration, 
for any initial guess. Therefore, in this case, solution of the linear system in the standard form 
(Eq. (21)) and solution in the incremental form (Eqs. (22) and (23)) are equivalent. 

More^ generally, however, the matrix [A] is not necessarily equal to the matrix [A], The 
matrix [A] can be any convenient approximation of the matrix [A] with the restriction that [A] 
must approximate [A] well enough so that the two-step iterative procedure (Eqs. (22) and (23)) 
converges (or, at the very least, can be forced to converge by including a strategy such as under- 
relaxation). Simply stated, [A] should capture the essence of [A]. Furthermore, because the 
equations have been cast in “delta” form, the incremental method produces the unique solution 
of Eq. (21), {Z* }, if convergent. In this formulation, the purpose of the left-hand side operator 
is to drive the right-hand side vector to zero. The final converged solution, {Z* }, depends only on 
the terms on the right-hand side of Eq. (22), and thus it is emphasized here that approximations 
to any of these terms, including the matrix [A], will produce erroneous final results. 
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In principle, the linear system of Eq. (22) can be solved either directly or iteratively, at each 
m th iteration level. If a direct method is chosen, only a single LU factorization of the coefficient 
matrix, [A], is needed, where the LU factorization is then reused for an unlimited number of 
iterations, including when multiple solutions of Eq. (21) are sought for different values of the 
vector, {B}. If the coefficient matrix, [A], is too large, an iterative algorithm will be the only 
recourse because of computer storage limitations. 

With the choice of an iterative algorithm, an “inner” iteration index, ‘i\ is established at 
Step 1 (Eq. (22)), and the iteration cycle over Steps 1 and 2,Jiaving index ‘m\ becomes the 
“outer” iteration loop. If the left-hand side coefficient matrix, [A], is diagonally dominant, then 
convergence of the iterative method of choice over the index, i , is assured for each and every 
linear sub-problem at Step 1. In addition, overall convergence^ the procedure over the outer 
index, ‘m’, is assured if, as discussed previously, the matrix [A] is an adequate approximation 
of the matrix [A], and furthermore, if each linear sub-problem at Step 1 is converged to a 
sufficiently close tolerance (whatever that tolerance may be). 

As a simple example of the preceding discussion.Jf a conventional relaxation algorithm (one 
of many possibilities) is selected, then the matrix, —[A], is divided into two parts, that is 

-[A] = [M] + [N] (24) 

The iterative incremental strategy becomes 

Step 1 [M] j m,i Az) = [A]{Z m } + {B) - [N]{”'' i " 1 Az} (2J) 

i = 1,2,3,..., (imax)" 1 

Step 2 {Z m+1 } = {Z m } + { m , (imax) m AZ } (26) 

m= 1,2,3,.... 

where in the above, (irnax) m is the number of inner or sub-iterations required to converge the m* 
linear sub-problem at Step 1 to the desired tolerance. The particular choice of the splitting in Eq. 
(24) of the matrix, [A], is made judiciously, such that Eq. (25) can be repeatedly solved very 
efficiently in terms of CPU time and computer storage. The most popular choices of the splitting 
in Eq. (24) result in either the Jacobi or the Gauss-Seidel algorithms of either the point or the 
line relaxation types. The use of the “delta” form line Gauss-Seidel algorithm with an “inner” 
and “outer” loop is investigated in Ref. 40 in the solution of the nonlinear 2D flow equations. 

2.5 Incremental Solution of the Equations of Sensitivity Analysis 

Application of the fundamental incremental formulation for linear equation solving, Eqs (22) 
and (23), to the linear system of Eq. (13) for computing aerodynamic sensitivity derivatives, 

gives 


12 




where the coefficient matrix 


dR 


approximates the matrix 



and will be discussed subse- 


quently, in greater detail. The vector, js m j, henceforth called the sensitivity residual 

vector, represents the total derivative of the discrete (flow analysis) residual vector, Eq. (12), 
with respect to /?k- From Eqs. (13) and (29), clearly the sensitivity residual vector must be 
driven to zero in order to find the solution, J, of Eq. (13), which is of course the objective 
of the incremental strategy of Eqs. (27) and (28). Approximations must not be made to any 
terms in the sensitivity residual vector, taking particular care that a consistent treatment of all 
boundary conditions is included here, if the converged solution is to yield the correct (i.e., the 
consistently discrete) sensitivity derivatives. The final solution at convergence depends only on 
the terms of this right-hand side vector. 


It is proposed that a first-order accurate upwind spatial discretization of the inviscid terms 
is a suitable selection for use in the coefficient matrix, ||^j , of Eq. (27), as an approximation 
here to the higher-order accurate upwind discretization of these terms. It is believed intuitively 
that this approximation will be a successful choice, noting that this selection is also a common 
approximation of convenience which is successfully used in the coefficient matrix of the implicit 
time integration formulation, Eq. (10). It is most significant to note that by design, in this choice, 
the block-diagonal dominance is now obtained and maintained in the left-hand side coefficient 
matrix (Ref. 29) of Eq. (27). 

In this preliminary study, the feasibility of using this first-order accurate upwind approximate 
treatment of the inviscid terms is investigated in the example problems. Of principal concern, of 
course, is whether or not this particular approximation yields a sufficiendy accurate representation 
of these terms so that a convergent method results. However, if the proposed methodology is 
successful, as it is found to be in the subsequent example problems, then the door has been opened 
for the possible future inclusion of numerous additional “approximations of convenience” in the 
left-hand side coefficient matrix. Of particular interest in future studies, of course, would be 
some of the same previously noted approximations commonly included in the coefficient matrix 
of the implicit formulation for time integration (Eq. (10)) of the flow equations. In other words, 
typical existing CFD flow solvers (i.e., those which employ iterative delta form implicit time 
integration methods) might be adapted directly for use in solving the linear sensitivity equations. 
The feasibility of this proposal is confirmed in the appendix, where sample results are presented 
using the well-known spatially-split approximate factorization (AF) (Ref. 34) algorithm. 
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In the present preliminary examination of the proposed methodology, each linear sub- 
problem (i.e., Eq. (27)) is solved by direct LU factorization (followed by forward and backward 
substitution) using a conventional vectorized banded matrix solver (Ref. 33) which takes 
advantage of the fact (in terms of computational work and storage) that outside of the bandwidth, 
all of the elements are zero. A single complete sensitivity analysis requires a single LU 
factorization of the coefficient matrix, which is repeatedly reused in the forward and backward 
substitutions at each iteration over Eqs. (27) and (28), and for all design variables of interest. 
Note that the direct solution of Eq. (27) now requires only one-half of the computer storage 
of that which is required in the direct solution of the equations in standard form, Eq. (13), 
since the bandwidth of the coefficient matrix is cut in half by the use of the first-order upwind 
approximation. In addition, less computational work is required in the LU factorization of this 
coefficient matrix, and in the forward and backward substitutions (although only a single back- 
solving procedure is required per design variable for a direct method applied to the standard 
formulation). 

The strategy proposed above is described as a combined iterative/direct solver method. It is 
felt however that the algorithm remains a direct solver method in its essential character, because 
despite the “factor of two” reduction in computer storage requirements, it remains infeasible 
for extension to practical 3D flow problems. However, the method will enable a significantly 
larger problem to be done in 2D. The present methodology will become purely iterative in 
character (and thus in principle extendable to 3D) when, as illustrated in Eqs. (25) and (26), 
an iterative method replaces the present direct solution of each linear sub-problem of Eq. (27). 
As an example given in the appendix, the AF algorithm is used to efficiently solve Eq. (27) 
approximately at each m^ iteration (without the use of sub-iterations), resulting in a convergent 
overall method. It is noted that convergence of jterative methods over each linear sub-problem 
(i.e, over each “inner loop”) is assured, since is block-diagonally dominant. 

Finally, it is noted here that if the adjoint variable formulation for computing the sensitivity 
derivatives is preferred, application of the incremental formulation to the linear system of Eq. 
(19) for computing the adjoint variable vector, { Aj } , yields 

r ~-, T 

Step 1 - { m AAj} = {V m (Aj 11 )} (30) 

Step 2 {Af+‘} = {Ai"}+rAA i } (J1) 

m = 1, 2, 3, 

where 

{v-W} = [fjj 

The vector, jv m j, known here as the adjoint variable residual vector, must be driven to 
zero in order to find the solution, {Aj }, of Eq. (19), which is the objective of the incremental 
strategy of Eqs. (30) and (31). 
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3.0 Computational Results 


Aerodynamic sensitivity derivatives computed using the incremental formulation, Eqs. (27) 
and (28), are presented for two laminar example problems, and are compared with the same 
results reported in Ref. 23 for the identical example problems. In Ref. 23, these same 
sensitivity derivatives were computed using direct solver based methods applied to the standard 
formulation of Eq. (13). It is significant to note that at the outset of this study all attempts 
to solve these sensitivity equations in standard form using a conventional line Gauss-Seidel 
iteration method (Refs. 29, 30) for these two example problems diverged, despite efforts to 
force convergence through the use of successive line under-relaxation. This failure is attributed 
to the ill-conditioning of the coefficient matrix. 

3.1 Internal Flow - Double- Throat Nozzle Problem 

The first example problem is that of an internal flow through a double-throat nozzle, where 
the flow is accelerated from a Mach number on the inflow boundary of about 0.10, to a Mach 
number which exceeds 2.80 at some places on the outflow boundary. The Reynolds number, 
REl, is 100, based on a reference length, L, of one-half the height of the nozzle at the smaller of 
the two throats. Figure 1 illustrates the geometry and computational grid which is used, and Fig. 
2 depicts the Mach contours of the steady-state solution; both of these figures are taken from Ref. 
23, where more complete information is given. Other studies have been conducted involving 
the numerical computation of flow through the geometry of this nozzle, and are documented in 
Refs. 41, 42, and 43. 

The geometric shape is defined parametrically using analytical expressions which define the 
boundaries (i.e., the walls) of the nozzle. Within these analytical expressions, ten geometric 
shape design variables are defined, and hence these ten parameters also define the vector, /?. 
These ten design variables, fi\ through /? io, the analytical functions which define this geometric 
shape, and also the treatment of the grid sensitivity vectors, through j^j, are fully 

explained in Refs. 23 and 43. 

In Ref. 23, the sensitivity derivatives were computed (with respect to (3\ through /3 10 ) of 
the force coefficients, C x and C y ; these force coefficients are the integrated (over the lower 
surface) pressure and skin frictions coefficients, which have been resolved in the V and ’y’ 
directions, respectively. In this earlier study, these terms were calculated by direct solution of the 
aerodynamic sensitivity equations in standard form (Eq. (13)), where a single LU factorization 
was used in the back-solving operations for all ten design variables. Additionally in the previous 
work, the accuracy of the calculations was successfully validated using the method of “brute 
force” finite differences, and thus this consistency check is not repeated here. 

In Table 1, the sensitivity derivatives of C x and C y with respect to the first geometric shape 
design variable, are presented. The computed values of and } are presented 

here from the solution of the aerodynamic sensitivity equations in incremental form, where 
results are given for successively larger reductions in the average global error. Specifically, the 
sensitivity derivatives computed using the incremental method are given for a zero, one, two, 
three, and four orders-of-magnitude (OM) reduction in the L 2 norm of the sensitivity residual 
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vector, Eq. (29), which is also the right-hand-side of Eq. (27). In addition, the number of 
iterations (over the two-step scheme, Eqs. (27) and (28)) which were required to achieve each 
of these successive levels of convergence is also included in the table. In the last row of the 
table, the results which were obtained by direct solution of the standard form of the equations, 
taken from Ref. 23, also are given. Tables 2 through 10 show results similar to those shown 
in Table 1, except that sensitivity derivatives of C x and C y with respect to /?2 through /3j 0 , 


respectively, are presented. 

For this first example problem, from the results presented in these tables, it is verified that the 
diagonally dominant first-order accurate upwind spatial discretization of the inviscid terms in the 


that the 


matrix, j^j , of Eq. (27) is a sufficiently accurate approximation of the matrix, 
iterative incremental formulation for solving these equations is convergent It is noted that these 


dR 

(Ki 


results were obtained without the use of under-relaxation or any scheme to “force” the method to 
converge. The solutions appear to be fairly well converged after only a two OM reduction of the 
error, and the first four digits (at least) of these sensitivity derivatives do not change as the error is 
reduced from three to four OM. Most importantly, the expected result is noted (as a consistency 
check), that the “tightly” converged results obtained using the incremental formulation agree 
with the results of Ref. 23 which were obtained using the standard formulation. 


Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

dfli 

dC y 


0 OM* 

1 

-3.877 E+01 

-3.211 E+02 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

13 

-4.934 E+01 

-3.024 E+02 

2 OM 

20 

-4.925 E+01 

-3.024 E+02 

3 OM 

27 

-4.925 E+01 

-3.024 E+02 


4 OM 

33 

-4.925 E+01 

-3.024 E+02 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

-4.925 E+01 

-3.024 E+02 


Table 1 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, First Design Variable, 

*OM Refers to the number of Orders-of-Magnitude reduction in the average global error. 
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Strategy 

Error 

Number of 


dCy 

d0 2 

Used 

Reduction 

Iterations 

X 

h 


0 OM 

1 

-4.644 E+02 

+2.152 E+01 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

8 

-4.614 E+02 

+1.733 E+01 

2 OM 

15 

-4.614 E+02 

+1.742 E+01 


3 OM 

22 

-4.614 E+02 

+1.741 E+01 


4 OM 

33 

-4.6 1' E+02 

+1.741 E+01 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

-4.614 E+02 

+1.741 E+01 


Table 2 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Second Design Variable, 02 


Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d/? 3 

d Cy 

d& 


0 OM 

1 

+2.343E +02 

-3.655 E+01 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

11 

+2.284 E+02 

-2.616 E+01 

2 OM 

18 

+2.284 E+02 

-2.625 E+01 


3 OM 

24 

+2.284 E+02 

-2.625 E+01 


4 OM 

31 

+2.284 E+02 

-2.625 E+01 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

+2.284 E+02 

-2.625 E+01 


Table 3 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Third Design Variable, 03 
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Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d/? 4 

dCy 

dft 

— 


0 OM 

1 


+2.213 E+03 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

10 

-2.665 E+04 

+1.659 E+03 

2 OM 

17 

-2.665 E+04 

+1.665 E+03 

3 OM 

23 

-2.665 E+04 

+1.664 E+03 


4 OM 

31 

-2.665 E+04 

+1.664 E+03 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

-2.665 E+04 

+1.664 E+03 


Table 4 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Fourth Design Variable, @4 


Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d/? 5 

dCy 

d& 


0 OM 

1 

-8.334 E+Ol 

+7.905 E-01 


1 OM 

3 

-8.326 E+01 

+4.500 E-01 

Incremental Method, 
Eqs. (27), (28), (29) 

2 OM 

6 ^ 

-8.327 E+01 

+4.344 E-01 

3 OM 

26 

-8.327 E+01 

+4.370 E-01 


4 OM 

46 

-8.327 E+01 

+4.365 E-01 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

-8.327 E+01 

+4.365 E-01 


Table 5 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Fifth Design Variable, fis 
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Table 6 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Sixth Design Variable, 06 


Strategy 

Error 

Number of 

dC x 

dC y 

Used 

Reduction 

Iterations 




d/3 7 

d/? 7 


0 OM 

1 

+2.120 E+00 

-5.216 E-01 

Incremental Method, 
Eqs. (27), (28), (29) 

OM 

12 

+1.444 E+00 

-5.873 E+00 

2 OM 

18 

+1.414 E+00 

-5.877 E+00 


3 OM 

25 

+1.415 E+00 

-5.879 E+00 


4 OM 

31 

+1.415 E+00 

-5.879 E+00 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

+1.415 E+00 

-5.879 E+00 


Table 7 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Seventh Design Variable, 0j 
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Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d(3 8 

dCy 

d/3 8 


0 OM 

1 

-3.415 E-01 

+2.353 E+02 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

15 

+6.281 E+00 

+2.331 E+02 

2 OM 

21 

+6.236 E+00 

+2.331 E+02 


3 OM 

28 

+6.236 E+00 

+2.331 E+02 


4 OM 

33 

+6.236 E+00 

+2.331 E+02 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

+6.236 E+00 

+2.331 E+02 


Table 8 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Eighth Design Variable, /3g 


Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d/? 9 

dCy 

d/? 9 


0 OM 

1 

-1.366 E+00 

-2.382 E+01 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

12 

-2.153 E+00 

-2.082 E+01 

2 OM 

18 

-2.107 E+00 

-2.082 E+01 


3 OM 

25 

-2.107 E+00 

-2.081 E+01 


4 OM 

31 


-2.081 E+01 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

N/A 

-2.107 E+00 

-2.081 E+01 


Table 9 - Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, Ninth Design Variable, 0g 
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Strategy 

Used 

Error 

Reduction 

Number of 
Iterations 

dC x 

d/?io 

dCy 

d/?io 


0 OM 

1 

+9.750 E-02 

+1.144 E+01 

Incremental Method, 
Eqs. (27), (28), (29) 

1 OM 

7 

+3.988 E-01 

+1.157 E+01 

2 OM 

13 

+3.903 E-01 

+1.158 E+01 


3 OM 

20 

+3.886 E-01 

+1.158 E+01 


'4 OM 

26 

+3.886 E-01 

+1.158 E+01 

Standard Form, Direct 
Solution of f- ... (13) 

N/A 

N/A 

+3.886 E-01 

+1.158 E+01 


Table 10 - C omparison of Sensitivity Derivatives, Incremental 
and Standard Methods, Tenth Design Variable, /? 10 


Table 1 1 shows a comparison of total CPU times, where naturally the computational cost of 
the incremental method depends heavily on the “strictness” of the desired convergence tolerance. 
For only a two OM error reduction, the computational cost of the incremental and the standard 
formulations are approximately equal. However, a tightly converged (four OM error reduction) 
solution results in a factor of almost two greater computational cost for the incremental method 
in the present example problem. 


Strategy 

Error 

CPU Time 

Used 

Reduction 

(Seconds)* 


0 OM 

27 

Incremental Method, 


51 

2 OM 

68 

Eqs. (27), (28), (29) 


3 OM 

90 


4 OM 

113 

Standard Form, Direct 
Solution of Eq. (13) 

N/A 

66 


Table 1 1 - Comparison of Total CPU Times, Incremental and Standard Methods 
*A11 Calculations Performed on a Cray-2 Computer. 
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It is noted that the sensitivity derivatives presented in the first row of Tables 1 through 
10 (i.e., for a zero OM error reduction, which is one iteration of the incremental method) are 
exactly the values which would be computed by direct solution of the standard formulation, if 
the left-hand side coefficient matrix, §j| , of Eq. (13) were approximated using the matrix, 

||| . By comparison of these calculations in the first row with those in the last row of the 
tables (i.e., the actual results of the standard formulation), the significant error is seen which 
would be generated in the sensitivity derivatives if approximations of convenience such as this 
were introduced into the standard formulation of the equations. 

3.2 External Flow - NACA 4-Digit Airfoil Problem 

The second problem considered here is that of external flow over an isolated airfoil, and is 
identical to the second example problem of Ref. 23. There pertinent details are found, including 
the grid and boundary conditions used, as well as an explanation of the special treatment for the 
grid sensitivity terms. The numerical solution of this laminar flow problem is for a freestream 
Mach number. Moo = 0.70, Reynolds number, REl = 5000, and angle of attack a = 0.0°. The 
airfoil shape is the NACA 2412, where the profile is defined by polynomial expressions in terms 
of three parameters, which are maximum thickness, T=0.12, maximum camber, C=0.02, and 
location of maximum camber, L=0.40. These three parameters are defined here to be the design 
variables, and hence define the elements of the vector, $. 

In Ref. 23, sensitivity derivatives were computed (with respect to T, C, and L) for the 
lift (Cl) and drag (Cd) coefficients. These terms were calculated in this earlier work using a 
hybrid direct solver/conventional iterative approach in the solution of the sensitivity equations 
in standard form (Eq. (13)). That is, a single direct LU factorization was applied to the central 
bandwidth of the coefficient matrix; the relatively small number of implicit terms which fall 
outside this main bandwidth (some at extreme distances because of the “periodic” boundary 
conditions at the “wake-cut” of the C-mesh) were treated “explicitly,” i.e., on the right-hand side 
of the equations. Thus a conventional Richardson iterative cycle was established to account for 
the periodic boundary conditions. However, despite the relatively small number of terms which 
were treated explicitly, it was reported that because of the required use of the poorly conditioned 
higher-order accurate coefficient matrix, the iterative strategy was at first divergent, and the use 
of under-relaxation was necessary to force the procedure to converge. As in the first example 
problem, the accuracy of the final results was successfully verified in this earlier work by finite 
differences, and thus this consistency check is omitted here. 

In the present application of the incremental strategy to this identical airfoil problem, the 
elements falling outside the central bandwidth of the left-hand coefficient side matrix, , were 
simply neglected entirely. This of course constitutes the inclusion of a second approximation 
of convenience in this matrix, in addition to the first-order accurate upwind treatment of the 
inviscid terms. The analogous (but not identical) terms resulting from the C-mesh type periodic 
boundary conditions in the matrix, , are not and must not be neglected on the right-hand 
side of Eq. (27), if the final sensitivity derivatives are to be correct. However, the treatment 
of these periodic terms is explicit and straightforward since they are on the right-hand side of 
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the equations. The resulting incremental strategy is again found to be convergent in the present 
example problem, with it the need for under-relaxation or any scheme to force the convergence. 
As in the first example 4 :oblem, the method is implemented by a single direct LU factorization 

of the approximate coefficient matrix, |q- , which is repeatedly reused in all subsequent back- 
solving operations, for all iterations and design variables. 

Table 12 shows the computed sensitivity derivatives of Cl and Cd with respect to f3\ = 
T, >r successively larger reductions in the error, where the results of the present incremental 
formulation are compared with the results for the standard formulation, taken from Ref. 23. 
Tables 13 and 14 provide similar results except that derivatives with respect to /? 2 = C and /? 3 
= L, respectively, are computed. Note that in these tables, the convergence of each method 
is fairly good after a two OM reduction in the error, and excellent after three or four OM. 
In addition, the converged results of the standard and incremental formulations are seen to 
consistently agree with one another, as expected. Table 15 presents the number of iterations 
required to achieve each level of error reduction, for each design variable, where the incremental 
and standard formulations are compared. Finally, Table 16 compares the total CPU times which 
were required in the calculations using the incremental and standard forms. For the present 
problem, the incremental method is seen to be more efficient. 


Error 

Reduction 

Lift Sensitivity 

dC L dC L 
dft dT 

Drag Sensitivity 

dC D dC D 
dft dT 

Standard 

Incremental 

Standard 

Incremental 

OOM 

-9.334 E-01 

-2.467 E-01 

+4.723 E-01 

+1.226 E+00 

1 OM 

-2.589 E+00 

-2.939 E+00 

+4.267 E-01 

+4.353 E-01 

2 OM 

-3.117 E+00 

-3.126 E+00 

+3.972 E-01 

+3.938 E-01 

3 OM 

-3.126 E+00 

-3.126 E+00 

+3.939 E-01 

+3.938 E-01 

4 OM 

-3.126 E+00 

-3.126 E+00 

+3.938 E-01 

+3.938 E-01 


Table 12 — Comparison of Sensitivity Derivatives, Incremental and 
Standard Methods, /?i = T (Maximum Thickness) 
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Error 

Reduction 

Lift Sensitivity 

dC L dC L 
d/? 2 dC 

Drag Sensitivity 

dC D dC D 
dfo dC 


Standard 

Incremental 

Standard 

Incremental 

OOM 

+5.206 E+00 

+4.706 E+00 

+3.429 E-01 

+6.778 E-01 

1 OM 

+4.175 E+00 

+2.973 E+00 

+3.780 E-01 

+3.785 E-01 

2 OM 

+3.988 E+00 

+3.976 E+00 

+3.663 E-01 

+3.640 E-01 

3 OM 

+3.968 E+00 

+3.968 E+00 

+3.603 E-01 

+3.603 E-01 

4 OM 

+3.968 E+00 

+3.968 E+00 

+3.603 E-01 

+3.603 E-01 


Table 13 — Comparison of Sensitivity Derivatives, Incremental 
and Standard Methods, / 3 2 = C (Maximum Camber) 



Standard Methods, = L (Location Of Maximum Camber) 
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Error 

Reduction 

Number of Iterations 
Standard 

Number of Iterations 
Incremental 

T 

C 

L 

T 

C 

L 

OOM 

1 

1 

1 

1 

1 

1 

1 OM 

13 

14 

5 

12 

14 

9 

2 OM 

64 

39 

24 

136 

45 

32 

3 OM 

219 

188 

49 

239 

203 

91 

4 OM 

300 

276 

195 

297 

269 

225 


Table 15 - Number of Iterations Required, Incremental and Standard Methods 


Error 

Reduction 

Total CPU Time (Seconds) 

Standard 

Incremental 

0 OM 

27 

10 

1 OM 

33 

15 

2 OM 

54 

41 

3 OM 

124 

89 

4 OM 

191 

127 


Table 16 - A Comparison of Total CPU Times, Incremental and Standard Methods 
4.0 Summary and Conclusions 

It has been shown herein that for the future development and application of efficient iterative 
methods for solving the aerodynamic sensitivity equations, there are significant advantages which 
can be exploited within the incremental formulation which are not seen in the standard form 
of these equations. These benefits are derived from the flexibility of the “delta” formulation, 
which allows any convenient approximation to be introduced into the left-hand side coefficient 
matrix (which operates on the “delta terms”) without affecting the final computed values of 
the sensitivity derivatives, provided the resulting sequence of successive iterations which are 
generated converges. Future efforts in algorithm development can now be directed at solving the 
sensitivity equations in delta form using conventional iterative strategies which are commonly 
employed in solving the nonlinear flow equations. The goal is to adapt existing CFD flow solvers 
in 2D and 3D with few or no changes to also solve the equations of aerodynamic sensitivity 
analysis. In this regard, preliminary results obtained to date are encouraging; in the appendix the 
feasibility of this proposal is confirmed in the example problems using a fully iterative solution 
process. 
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6.0 Appendix - Future Work, An Approximately Factored Method 

Having developed and successfully demonstrated an incremental formulation which is flex- 
ible in character for solving the sensitivity equations, future work in algorithm development for 
the iterative solution of these equations will seek to adapt iterative strategies which are com- 
monly used in the implicit time integration of the flow equations. To this end, a false time term, 
which is the diagonal matrix, * * s added to the left-hand side coefficient matrix, , of 
Eq. (27). This “time” term diagonal matrix is of course found in the implicit time integration 
formulation of Eq. (10). 

The addition of this false time term to each element on the diagonal of the coefficient matrix 
in Eq. (27) is equivalent to the use of under-relaxation in the two-step incremental formulation 
of Eqs. (27) and (28). Then, for small to moderate time steps, the resulting linear system of 
Eq. (27) may be very efficiently solved (approximately) at each iteration (i.e., at each false time 
step) using the spatially-split approximate factorization algorithm (AF) of Ref. 34. This basic 
algorithm, which has many variations, is well known as a common strategy found in 2D and 3D 
CFD flow codes for the efficient approximate solution of Eq. (10) at each time step. 

With the introduction of the false time term to the elements on the diagonal, and the resulting 
factorization error which is associated with the AF algorithm at each iteration (all in addition 
to the error of the approximate first-order accurate upwind treatment of the inviscid terms), it 
was not known a priori whether the resulting approximate coefficient matrix operator on the 
left-hand side of Eq. (27) would be a convergent method for solving these equations. However, 
the proposed AF strategy has been found to be convergent in application to the two previously 
explained example problems of this study. That is, the algorithm was successfully used to 
produce a four OM reduction in the average error (as defined previously) for the double-throat 
nozzle problem and the airfoil problem. 

Using a constant Courant number of 10 for each cell in the computational grid (i.e., using 
local false “time-stepping”), for the double-throat nozzle example. Table 17 shows the computed 
sensitivity derivatives of C x and C y (along the lower wall) with respect to /?, through /? 10 , 
following the four OM reduction in the average error, where the number of iterations required 
by this algorithm to achieve this level of convergence is reported for each design variable. As 
expected, these results shown here agree very well with those reported earlier for this example 
problem, except some of the very small sensitivity derivatives show minor discrepancies which 
prove to disappear when the AF method is used to reduce the error to a stricter tolerance than 
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the four OM shown here. Table 18 presents a comparison of the total CPU time required in this 
example using the AF method compared to the CPU times shown earlier for the other methods. 


Design 

Variable 

Number of 
Iterations 

Sensitivity 
of C x 

Sensitivity 

of Cy 

01 

335 

-4.926 E+01 

-3.024 E+02 

02 

277 

-4.614 E+02 

+1.741 E+01 

03 

242 

+2.284 E+02 

-2.625 E+01 

04 

276 

-2.665 E+04 

+1.664 E+03 

05 

259 

-8.327 E+01 

+4.370 E-01 

06 

278 

-1.778 E-02 

+1.428 E+02 

07 

225 

+1.414 E+00 

-5.881 E+00 

08 

317 

+6.233 E-tOO 

+2.331 E+02 

09 

280 

! 

-2.109 E+00 

-2.081 E+01 

010 

243 

+3.881 E-01 

+1.158 E+01 


Table 17 - Double-Throat Nozzle Problem, Approximately Factored 
(AF) Incremental Method, Four OM Error Reduction 


Strategy 

Used 

Total CPU Time 
(Seconds) 

Incremental, AF 
Solver, (4 OM) 

144 

Incremental, Direct 
Solver (4 OM) 

113 

Standard Form, 
Direct Solution 

66 


Table 18 - Double-Throat Nozzle Problem, Comparison of Total CPU Times 
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Using a constant Courant number for each cell of 20 for the airfoil problem. Table 19 shows 
the computed sensitivity derivatives of Cl and Co with respect to T, C, and L, and the number 
of iterations required by the AF method are also given. As expected, the computed sensitivity 
derivatives here are in excellent agreement with those reported previously for this problem. Table 
20 is a summary of the total CPU times required in this example, comparing the present method 
with the previously presented results. 


Design 

Variable 

Number of 
Iterations 

dC L 

dC D 

dT, C, L 

dT, C, L 

T 

466 

-3.126 E+00 

+3.938 E-01 

C 

428 

+3.968 E+00 

+3.603 E-01 

L 

360 

-1.816 E-02 

-3.290 E-03 


Table 19 - NACA 2412 Airfoil Problem, Approximately Factored 
(AF) Incremental Method, Four OM Error Reduction 


Strategy 

Total CPU Time 

Used 

(Seconds) 

Incremental, AF 

50 

Solver, (4 OM) 


Incremental, Direct 

127 

Solver (4 OM) 


Standard Form, Hybrid 

191 

Direct/Iterative (4 OM) 



Table 20 - NACA 2412 Airfoil Problem, Comparison of Total CPU Times 

The preceding results are encouraging, and demonstrate the feasibility of the proposed 
methods. Much work remains in selecting and refining the most efficient algorithms and 
convergence accelerations methods (such as multigrid, for example) for use in the solution 
of the aerodynamic sensitivity equations in 2D and 3D. 
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